package scu.maqiang.matlabcases;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class LShapeEigenvalue {
	public static void main(String[] args) {
		String meshFile = "matlabMeshes\\Lshape.msh";
		Mesh2T3 mesh2D = new Mesh2T3().readMesh(meshFile, MeshFormat.msh);
		FES2T31 fs = new FES2T31(mesh2D);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		SRMatrix M = new SRMatrix(fs.GetNdof());
		fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
		fs.assembleMass(new double[]{1.0}, BVPType.CONSIST_MASS, M);
		double[] x = MVO.rand(fs.GetNdof());
		EigenSSolver solver = new EigenSSolver();
		int n = 20;
		double[] eigenvalues = new double[n];
		double[][] eigenvectors = new double[n][fs.GetNdof()];
		fs.applyBCZero(A, x, Direct.X, BCType.RCE, 1, 2, 3, 4, 5, 6);
		solver.inversePowerTran_DSS(A, M, x, 2.0, n, eigenvalues, eigenvectors);
		System.out.println("Eigenvalues: ");
		System.out.println(MVO.toString(eigenvalues));
		mesh2D.toTecplot("matlabResults\\LShapeEigenvectors.dat", eigenvectors);
	}

}
